Autocorrelation

Possible fixes (that don’t require throwing out data )

“Random” vs. Fixed

(the terminology here is inconsistent and contested, but no one has come up with better words yet)

• Fixed Effects Models give every group its own y-intercept. The effects of “Texas” vs. “Tennessee” are assumed to be fixed and unrelated.

  • In this sense, you’re usually assuming a fixed effect any time you run a standard OLS model.

• Random Effects Models assume that group intercepts are random draws from a common distribution. They can vary just like any other random variable, but they vary in a predictable way

An example: students in classrooms

flowchart LR
Z[School] --> A
Z --> G
Z --> L
A[Classroom 1] --> B(Student A)
A --> C(Student B)
A --> E(Student D)
G[Classroom 2] --> H(Student H)
G --> I(Student I)
G --> K(Student K)
L[Classroom 3] --> M(Student M)

An example: students in classrooms

You calculate some averages and your results look like this:

An example: students in classrooms

So we could assume something like this:

\[ SAT = \text{teacher effect} + b\_\text{studying effect} \times \text{study time} + \epsilon \] \[ \epsilon \sim \mathcal{N}(0, \sigma) \]

But we could assume a model with some kind of a “teacher effect” that also had a normal distribution like the error term:

\[ SAT = (\text{teacher effect}) + b\_\text{studying effect} \times \text{study time} + \epsilon\]

\[ \text{teacher effect} \sim \mathcal{N}(0, \tau) \]

Random Effects: pooling

One way to think about a random effect is that it represents a compromise between “no pooling” and “complete pooling” of regression results.

Random Effects: pooling

A partial pooling estimate favors “no pooling” estimate if group \(j\) is large and low-variance, and favors the “complete pooling” case if group \(j\) is small or high variance.

\[ \text{Estimate of random intercept } \alpha_j \approx \frac{\frac{n_j}{\sigma^2_j}}{\frac{n_j}{\sigma^2_j} + \frac{1}{\sigma^2_j}} (\bar{y_j} - \beta\bar{x}_j) + \frac{\frac{1}{\sigma^2_a}}{\frac{n_j}{\sigma^2_y} + \frac{1}{\sigma^2_a}} \mu_a \]

Random Effects: pooling

Random Effects: pooling

• Random effects models will often give more plausible estimates for cases with small sample sizes - because they “pool” extreme values with small samples.

  • Think of it like estimating a batting average for a player with 2 at-bats. Is plausible they’re batting 100? Or should we assume they’re around the league average until we get more information?

Random Effects

library(fixest)
# county level election results

fmodel<-feols(perc_gop ~ Income + White + AfAm + Hisp  | state_name, data=counties)
regmodel<-feols(perc_gop ~ Income + White + AfAm + Hisp , 
                data =counties, cluster ~ state_name)

mlist<-list("No fixed effects" = regmodel, 
            'Fixed effects' = fmodel)
modelsummary(mlist, 
             estimate  = "{estimate}",  
             statistic = c("conf.int", 'std.error'),
             conf_level = .95,
             coef_rename = coefnames,
             coef_omit ='Intercept',
             gof_map = c('nobs','aic','bic', 'vcov.type',
                         'FE: state_name'
             )
)

Random Effects

To specify a random effect, we can use the lme4 package. The (1|state_name) says that I want to estimate random intercepts for each state:

library(lme4)
# fixed effects version for comparison:
fixefmodel<-fixest::feols(perc_gop ~ Income + White + AfAm + Hisp  | state_name, 
                          data=counties)

ranefmodel<-lmer(perc_gop ~ Income + White + AfAm + Hisp + (1|state_name),
               data=counties
               )

Random Effects

Random Effects: additional covariates

With fixed effects models, we can’t include any controls for group-level covariates because they’re already in the “fixed” effect:

library(rvest)
page<-read_html('https://en.wikipedia.org/wiki/List_of_U.S._states_by_date_of_admission_to_the_Union')

join_date<-page|>
  html_element(css='table')|>
  html_table()|>
  select(2:4)|>
  mutate(datestring = str_extract(
    `Date(admitted or ratified)`, "([A-Z][a-z]+ [0-9]{1,2}, [0-9]{4})"),
         date = mdy(datestring)
  )
counties<-counties|>
  left_join(join_date, by=join_by(state_name==State))|>
  mutate(join_year = year(date) - 1787
         
         )|>
  drop_na(join_year)

fixefmodel<-fixest::feols(perc_gop ~ Income + White + AfAm + Hisp + join_year  | state_name, 
                          data=counties)

fixefmodel

OLS estimation, Dep. Var.: perc_gop
Observations: 3,114
Fixed-effects: state_name: 50
Standard-errors: Clustered (state_name) 
        Estimate Std. Error   t value   Pr(>|t|)    
Income -0.257136   0.025502 -10.08291 1.5294e-13 ***
White   0.591932   0.033514  17.66219  < 2.2e-16 ***
AfAm   -0.248722   0.056598  -4.39450 5.9451e-05 ***
Hisp    0.274869   0.061121   4.49710 4.2332e-05 ***
... 1 variable was removed because of collinearity (join_year)
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 8.35604     Adj. R2: 0.713165
                Within R2: 0.567463

Random Effects: additional covariates

But random effects models can still estimate effects for group-level characteristics:

ranefmodel<-lmer(perc_gop ~ Income + White + AfAm + Hisp + join_year + (1|state_name),
               data=counties
               )

summary(ranefmodel)

Linear mixed model fit by REML ['lmerMod']
Formula: perc_gop ~ Income + White + AfAm + Hisp + join_year + (1 | state_name)
   Data: counties

REML criterion at convergence: 22348.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.1757 -0.4745  0.0931  0.5828  3.3831 

Random effects:
 Groups     Name        Variance Std.Dev.
 state_name (Intercept) 113.01   10.63   
 Residual                71.06    8.43   
Number of obs: 3114, groups:  state_name, 50

Fixed effects:
            Estimate Std. Error t value
(Intercept) 31.27089    3.36819   9.284
Income      -0.25919    0.01074 -24.124
White        0.59241    0.02327  25.463
AfAm        -0.24204    0.02853  -8.482
Hisp         0.27559    0.02739  10.061
join_year    0.06814    0.03294   2.068

Correlation of Fixed Effects:
          (Intr) Income White  AfAm   Hisp  
Income    -0.276                            
White     -0.693  0.072                     
AfAm      -0.630  0.154  0.830              
Hisp      -0.584  0.061  0.816  0.691       
join_year -0.576  0.018  0.095  0.105  0.058

Random Effects: random slopes

We can also estimate random slopes using this same basic approach.

ranefmodel<-lmer(perc_gop ~  White + AfAm + Hisp + join_year + (Income |state_name),
               data=counties
               )

library(marginaleffects)
nd<-datagrid(model=ranefmodel, 
             state_name = c("Alabama", "Wyoming", "Maryland", "Virginia", "Nevada", "North Carolina"), 
             "Income" = seq(17, 170, by=20))

plot_predictions(ranefmodel, newdata=nd, by=c('Income', 'state_name')) +
  theme_bw() +
  labs(color ='State', fill='State',  y='% GOP vote')

What to use

• Fixed Effects model are best when:
  • You want to control for un-modeled differences across groups
  • You think groups are fundamentally “different”
  • You don’t want to control for any group-level characteristics
• Random Effects models are preferable when:
  • There are lots of small groups with very few members
  • You’re less worried about unmeasured heterogeneity between clusters
  • You want to include characteristics across different “levels” or even allow the slopes to vary or certain observations

Time

• OLS assumes that data are independent, but your height today is not independent of your height yesterday. This violates the assumption of independence of errors.
• Moreover, time be a confounding variable: if two things are impacted by time, they’ll appear correlated in a regression that fails to account for a trend.

Potential fixes

• Correcting the standard errors uses the similar methods to the ones we’ve discussed elsewhere.

• We could control for time linearly or with a polynomial or spline function.

• We could treat time as a factor and include it as a fixed effect.

• Would could treat time as a random effect term.

Potential fixes: time coefficients

model1<-lm(v2x_polyarchy ~ log(e_pop), data=vdem)

model2<-lm(v2x_polyarchy ~ log(e_pop) + year, data=vdem)

model3<-lm(v2x_polyarchy ~ log(e_pop) + splines::ns(year), data=vdem)


mlist<-list(model1, model2, model3)

modelsummary(mlist)

Potential fixes: fixed effects

If we have panel data (multiple groups over time) we can include year as a fixed effect. This has the effect of controlling for anything that varies for all groups over time.

fmodel1<-feols(v2x_polyarchy ~ log(e_pop) | year, data=vdem)
modelsummary(fmodel1)

Potential fixes: fixed effects

Since we have both time periods AND groups, we can also add country fixed effects:

fmodel1<-feols(v2x_polyarchy ~ log(e_pop) | year, data=vdem)
fmodel2<-feols(v2x_polyarchy ~ log(e_pop) | country_name + year, data=vdem)

mlist<-list(fmodel1, fmodel2)

modelsummary(mlist)

Potential fixes: differencing

• If data depends only on t-1, taking the first difference should address it

• difference = current_value - prior_value

Potential fixes: lagged variable model

Another option is to include a lagged DV variable itself as a covariate. This is typically referred to as an autoregressive model.

ar_model<-vdem|>
  filter(country_name=='Paraguay')|>
  arrange(year)|>
  mutate(lag_polyarchy = dplyr::lag(v2x_polyarchy),
         lag_gdp = dplyr::lag(e_gdppc)
         
         )|>
  filter(year>=1995)|>
  lm(v2x_polyarchy ~ lag_polyarchy  + log(e_gdppc)  + log(lag_gdp) , data=_)

Other fixes

• Including a fixed effect for time can work if you have multiple observations for different periods or if you only suspect a cyclical trend

• More complicated autocorrelation structures can be modeled with a regression, and then OLS can be used on the residuals from that model